function [drIonoG, TIono] = ionogrpdelay_klob(sysTime, ionoPar, phiU, lambdaU, E, A, kf, sysType, c)
% IONOGRPDELAY_KLOB 用Klobuchar模型计算单频电离层群延迟
%
% Tips
% # 本函数计算的电离层群延迟适用于伪码测距，对于载波相位测距的电离层超前需要取负号
%
% Input Arguments
% # sysTime: 标量，用户接收机计算出来的导航系统时（系统的周内秒），值域：0~604800，单位：s
% # ionoPar: 电离层延迟改正模型参数，结构体，包含以下域：
%     alpha: 长度为4的向量，分别为导航电文中的α0~α3参数，单位分别为s、s/π、s/π^2、s/π^3
%     beta: 长度为4的向量，分别为导航电文中的β0~β3参数，单位分别为s、s/π、s/π^2、s/π^3
% # phiU: 标量，用户接收机的大地纬度φu，值域：-π/2~π/2，单位：rad
% # lambdaU: 标量，用户接收机的大地经度λu，值域：-π~π，单位：rad
% # E: 标量，用户接收机到导航卫星之间的高度角，值域：0~π/2，单位：rad
% # A: 标量，用户接收机到导航卫星之间的方位角，值域：0~2π，单位：rad
% # kf: 标量，频率因子，无量纲
% # sysType: GNSSSysType枚举类型，GNSS系统类型
% # c: 标量，光在真空中的传播速度，单位：m/s
%
% Output Arguments
% # drIonoG: 标量，电离层群延迟对应的路径长度，单位m
% # TIono: 标量，电离层群延迟时间（针对GPS L1或北斗B1频点，未考虑频率因子），单位s
%
% Assumptions and Limitations
% # 本函数中使用的电离层修正算法是Klobucher模型，仅适用于地面用户使用
%
% References:
% # 《应用导航算法工程基础》“单频接收机的电离层延迟改正”

DC = 5.0e-9;
phi = 50400;

switch sysType
    case GNSSSysType.GPS
        %% GPS系统
        % 计算用户接收机至地心与电离层穿刺点到地心之间的夹角
        Psi = 0.0137/(E/pi+0.11) - 0.022;
        
        % 计算电离层穿刺点的地理纬度
        phiIInPi = phiU/pi + Psi*cos(A);
        if phiIInPi > 0.416
            phiIInPi = 0.416;
        elseif phiIInPi < -0.416
            phiIInPi = -0.416;
        end
        
        % 计算电离层穿刺点的地理经度
        lambdaIInPi = lambdaU/pi + Psi*sin(A)/cos(phiIInPi*pi);
        
        %计算电离层穿刺点的地磁纬度
        phiM = phiIInPi + 0.064*cos((lambdaIInPi-1.617)*pi);
        
        % 计算电离层延迟变化的最大振幅
        AMP = ionoPar.alpha(1) + phiM*(ionoPar.alpha(2) + phiM*(ionoPar.alpha(3) + phiM*ionoPar.alpha(4)));
        if AMP < 0
            AMP = 0;
        end
        A = AMP;
        
        % 计算电离层延迟变化的最大周期
        PER = ionoPar.beta(1) + phiM*(ionoPar.beta(2) + phiM*(ionoPar.beta(3) + phiM*ionoPar.beta(4)));
        if PER < 72000
            PER = 72000;
        end
        P = PER;
        
        % 计算倾斜因子
        F = 1.0 + 16*((0.53-E/pi)^3);
        
        % 计算穿刺点地方时
        t = 4.32e4*lambdaIInPi + sysTime;
        t = wrap(t, 0, 86400);
    case {GNSSSysType.GLONASS, GNSSSysType.Galileo}
        error('ionogrpdelay_klob:unsupportedgnsssys', 'GLONASS and Galileo systems are not supported.');
    case GNSSSysType.BeiDou
        %% 北斗系统
        RInKm = 6378;
        hInKm = 375;
        
        % 计算用户与穿刺点地心张角
        Psi = pi/2 - E - asin(RInKm*cos(E)/(RInKm+hInKm));
        
        %计算电离层穿刺点的地理纬度
        phiM = asin(sin(phiU)*cos(Psi) + cos(phiU)*sin(Psi)*cos(A));
        lambdaM = lambdaU + asin(sin(Psi)*sin(A)/cos(phiM));
        
        % 计算电离层延迟变化的最大振幅
        A2 = ionoPar.alpha(1) + abs(phiM)*(ionoPar.alpha(2) + abs(phiM)*(ionoPar.alpha(3) + abs(phiM)*ionoPar.alpha(4)));
        if A2 < 0
            A2 = 0;
        end
        A = A2;
        
        % 计算电离层延迟变化的最大周期
        A4 = ionoPar.beta(1) + abs(phiM)*(ionoPar.beta(2) + abs(phiM)*(ionoPar.beta(3) + abs(phiM)*ionoPar.beta(4)));
        if A4 >= 172800
            A4 = 172800;
        elseif A4 < 72000
            A4 = 72000;
        end
        P = A4;
        
        % 计算倾斜因子
        F = 1 / sqrt(1-(RInKm*cos(E)/(RInKm+hInKm))^2);
        
        % 计算穿刺点地方时（按GPS ICD文件计算，北斗ICD没有相关内容）
        t = 4.32e4*lambdaM/pi + sysTime;
        t = wrap(t, 0, 86400);
    otherwise
        error('ionogrpdelay_klob:unknowngnsssys', 'GNSS system is not GPS, GLONASS, Galileo or BeiDou.');
end

%% 计算电离层群延迟
X = 2*pi*(t-phi) / P;
if abs(X) < pi/2
    if sysType == GNSSSysType.GPS
        TIono = F * (DC+A*(1-(X^2)*(1/2-(X^2)/24)));
    else
        TIono = F * (DC+A*cos(X));
    end
else
    TIono = F * DC;
end
drIonoG = kf * c * TIono;
end